---
title: "pa2_analyses"
output: word_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r load libraries}
library(readr)
library(reshape2)
library(plyr)
library(ggplot2)
library(mediation)
library(matrixStats)
library(tidyverse)

```
###BASELINE ANALYSES####
```{r import}

#load cleaned csv from matlab
fulldata <- read_csv("~/Desktop/Studies/PunishmentAversion/Study2/baseline/pa2_baseline_fulldata_foranalyses.csv",col_names = TRUE)

```
#Clean baseline data
```{r clean_data}

#Remove rows with any missing data
cleandata <- fulldata[complete.cases(fulldata), ]

#Remove any subject who did not pass the qualitative comprehension
cleandata <- cleandata[which(cleandata$qual_comp_score == 1), ]

#Specify that the variables are factors
cleandata$condition<-factor(cleandata$condition,labels=c("action","outcome", "punishment"))

#separate out the scenarios
actiondata=cleandata[which(cleandata$condition=="action"),]
outcomedata=cleandata[which(cleandata$condition=="outcome"),]
punishmentdata=cleandata[which(cleandata$condition=="punishment"),]
```
#Calculate baseline averages
```{r calculate averages}
#get column numbers for first and last action columns in the dataset
first_col <- which(colnames(cleandata) == "shirt")
last_col <- which(colnames(cleandata) == "restrain_20")

#Calculate the average judgment for each punishment for each condition (saved as vectors)
actionaverage = colMeans(actiondata[first_col:last_col])

outcomeaverage = colMeans(outcomedata[first_col:last_col])

punishmentaverage = colMeans(punishmentdata[first_col:last_col])

#create averages dataframe (combine vectors into single data frame)
averages.data <- data.frame(actionaverage, outcomeaverage, punishmentaverage)

#create a variable that specifies the scenario
averages.data <- rownames_to_column(averages.data, var = "scenario")

```

#Baseline Plots
```{r basic_plots}

#plot action vs. outcome
ggplot(averages.data,aes(y=actionaverage,x=outcomeaverage))+
  geom_point()+
  geom_smooth(method=lm)


#plot punishment predicted by action
ggplot(averages.data,aes(y=punishmentaverage,x=actionaverage))+
  geom_point()+
  geom_smooth(method=lm)

#plot punishment predicted by outcome
ggplot(averages.data,aes(y=punishmentaverage,x=outcomeaverage))+
  geom_point()+
  geom_smooth(method=lm)

```
#Baseline Analyses

```{r analyses}

#test correlation between action and outcome
mosaic::cor.test(actionaverage~outcomeaverage, data = averages.data)

#predict punishment score from action and outcome scores
mod1 <- (lm(punishmentaverage ~ outcomeaverage + actionaverage))
summary(mod1)

#test the relative effect of outcome and average
library(car)
linearHypothesis(mod1, c("actionaverage = outcomeaverage"))


```

###INGROUP/OUTGROUP ANALYSES###

```{r import}

#load cleaned csv from matlab
ingroup_outgroupdata <- read_csv("~/Desktop/Studies/PunishmentAversion/Study2/ingroup_outgroup/pa2_ingroup_outgroup_fulldata_foranalyses.csv", col_names = TRUE)

```
#Clean ingrou/outgroup data
```{r clean_data}

#Remove rows with any missing data
clean_inout_data <- ingroup_outgroupdata[complete.cases(ingroup_outgroupdata), ]

#Remove any subject who did not pass the qualitative comprehension
clean_inout_data <- clean_inout_data[which(clean_inout_data$qual_comp_score == 1), ]

#Specify that the variables are factors
clean_inout_data$condition<-factor(clean_inout_data$condition,labels=c("ingroup","outgroup"))

#separate out the scenarios
ingroup_punishmentdata=clean_inout_data[which(clean_inout_data$condition=="ingroup"),]
outgroup_punishmentdata=clean_inout_data[which(clean_inout_data$condition=="outgroup"),]
```

#Calculate ingroup/outgroup averages
```{r calculate averages}
#get the column numbers for the first and last scenario columns
first_col <- which(colnames(clean_inout_data) == "shirt")
last_col <- which(colnames(clean_inout_data) == "restrain_20")

#Calculate the average judgment for each punishment for each condition (saved as vectors)
ingroup_punishmentaverage = colMeans(ingroup_punishmentdata[first_col:last_col])
outgroup_punishmentaverage = colMeans(outgroup_punishmentdata[first_col:last_col])

#create final averages dataframe make 2 and join them together
ingroupdata <- data.frame(actionaverage, outcomeaverage,ingroup_punishmentaverage)
colnames(ingroupdata)[3] <- c("punishmentaverage")
#set the condition to 1 when it is ingroup condition
ingroupdata$condition = 1

#create a variable that specifies the scenario
ingroupdata <- rownames_to_column(ingroupdata, var = "scenario")


outgroupdata <-data.frame(actionaverage, outcomeaverage,outgroup_punishmentaverage)
colnames(outgroupdata)[3] <- c("punishmentaverage")
#set the condition to -1 when it is outgroup
outgroupdata$condition = -1

#create a variable that specifies the scenario
outgroupdata <- rownames_to_column(outgroupdata, var = "scenario")

#combine into a full dataset
inout_averages.data <- rbind(ingroupdata,outgroupdata)

#mean center numerical variables
inout_averages.data$actionavg_cent<-inout_averages.data$actionaverage-(mean(inout_averages.data$actionaverage))
inout_averages.data$outcomeavg_cent<-inout_averages.data$outcomeaverage-(mean(inout_averages.data$outcomeaverage))
inout_averages.data$punishmentavg_cent<-inout_averages.data$punishmentaverage-(mean(inout_averages.data$punishmentaverage))

#create interaction variables
inout_averages.data$aa_x_cond <- as.numeric(inout_averages.data$condition)*inout_averages.data$actionavg_cent
inout_averages.data$oa_x_cond <- as.numeric(inout_averages.data$condition)*inout_averages.data$outcomeavg_cent

#specify that categorical variables are factors
inout_averages.data$scenario <- factor(inout_averages.data$scenario)
inout_averages.data$condition <- factor(inout_averages.data$condition, labels = c("outgroup", "ingroup"))

```
#Ingroup/outgroup descriptives
```{r descriptives}
punishment_desc <- ddply(inout_averages.data, c("condition"), summarise,
               N    = length(punishmentaverage),
               mean = mean(punishmentaverage),
               sd   = sd(punishmentaverage),
               se   = sd / sqrt(N)
)
print(punishment_desc)

ingroup_data<-inout_averages.data[which(inout_averages.data$condition=="ingroup"),]
outgroup_data<-inout_averages.data[which(inout_averages.data$condition=="outgroup"),]

mosaic::cor(punishmentavg_cent~actionavg_cent,data=ingroup_data)
mosaic::cor(punishmentavg_cent~actionavg_cent,data=outgroup_data)

mosaic::cor(punishmentavg_cent~outcomeavg_cent,data=ingroup_data)
mosaic::cor(punishmentavg_cent~outcomeavg_cent,data=outgroup_data)
```
#Ingroup/outgroup analyses
```{r linear models}

# full model
full_mod <- lme4::lmer(punishmentavg_cent ~ actionavg_cent + outcomeavg_cent + condition + aa_x_cond + oa_x_cond + (1|scenario),data=inout_averages.data)

summary(full_mod)

#to test the hypothesis that the actionavg beta is significantly larger than the outcomeavg beta
library(car)
linearHypothesis(full_mod, c("actionavg_cent = outcomeavg_cent"))

#model without action main effect
NoAction_mod <- lme4::lmer(punishmentavg_cent ~ outcomeavg_cent + condition + aa_x_cond + oa_x_cond + (1|scenario),data=inout_averages.data)

summary(NoAction_mod)
anova(full_mod,NoAction_mod)

#model without outcome main effect
NoOutcome_mod <- lme4::lmer(punishmentavg_cent ~ actionavg_cent + condition + aa_x_cond + oa_x_cond + (1|scenario),data=inout_averages.data)

summary(NoOutcome_mod)
anova(full_mod,NoOutcome_mod)

#model without condition main effect
NoCondition_mod <- lme4::lmer(punishmentavg_cent ~ actionavg_cent + outcomeavg_cent + aa_x_cond + oa_x_cond + (1|scenario),data=inout_averages.data)

summary(NoCondition_mod)
anova(full_mod,NoCondition_mod)

#model without action x condition interaction
NoActxCond_mod <- lme4::lmer(punishmentavg_cent ~ actionavg_cent + outcomeavg_cent + condition + oa_x_cond + (1|scenario),data=inout_averages.data)

summary(NoActxCond_mod)
anova(full_mod,NoActxCond_mod)

#model without outcome x condition interaction
NoOutxCond_mod <- lme4::lmer(punishmentavg_cent ~ actionavg_cent + outcomeavg_cent + condition + aa_x_cond + (1|scenario),data=inout_averages.data)

summary(NoOutxCond_mod)
anova(full_mod,NoOutxCond_mod)

```
#Ingroup/outgroup plots
```{r new graphs}

# action vs punishment + outcome vs punishment divided by condition

# punishment vs action 
ggplot(inout_averages.data,aes(y=punishmentaverage,x=actionaverage))+
  geom_point()+
  facet_wrap(~condition)+
  geom_smooth(method=lm)

# punishment vs outcome 
ggplot(inout_averages.data,aes(y=punishmentaverage,x=outcomeaverage))+
  geom_point()+
  facet_wrap(~condition)+
  geom_smooth(method=lm)

#action vs. outcome
ggplot(inout_averages.data,aes(y=actionaverage,x=outcomeaverage))+
  geom_point()+
  geom_smooth(method=lm)

#bar graphs

#punishment
ggplot(data=inout_averages.data, aes(x=scenario, y=punishmentaverage, fill=condition)) +
  geom_bar(stat="identity", width=.9, position=position_dodge()) + scale_fill_brewer(palette="Set2")


```
###SELF/OTHER ANALYSES###

```{r import}

#load cleaned csv from matlab
self_otherdata <- read_csv("~/Desktop/Studies/PunishmentAversion/Study2/self_other/pa2_self_other_fulldata_foranalyses.csv", col_names = TRUE)

```
#Clean self/other data
```{r clean_data}

#Remove rows with any missing data
clean_selfother_data <- self_otherdata[complete.cases(self_otherdata), ]

#Remove any subject who did not pass the qualitative comprehension
clean_selfother_data <- clean_selfother_data[which(clean_selfother_data$qual_comp_score == 1), ]

#Specify that the variables are factors
clean_selfother_data$condition<-factor(clean_selfother_data$condition,labels=c("self","other"))

#separate out the scenarios
self_punishmentdata=clean_selfother_data[which(clean_selfother_data$condition=="self"),]
other_punishmentdata=clean_selfother_data[which(clean_selfother_data$condition=="other"),]
```

#Calculate self/other averages
```{r calculate averages}
#get the column numbers for the first and last scenario columns
first_col <- which(colnames(clean_selfother_data) == "shirt")
last_col <- which(colnames(clean_selfother_data) == "restrain_20")

#Calculate the average judgment for each punishment for each condition (saved as vectors)
self_punishmentaverage = colMeans(self_punishmentdata[first_col:last_col])
other_punishmentaverage = colMeans(other_punishmentdata[first_col:last_col])

#create final averages dataframe make 2 and join them together
selfdata <- data.frame(actionaverage, outcomeaverage,self_punishmentaverage)
colnames(selfdata)[3] <- c("punishmentaverage")
#set the condition to 1 when it is self condition
selfdata$condition = 1

#create a variable that specifies the scenario
selfdata <- rownames_to_column(selfdata, var = "scenario")


otherdata <-data.frame(actionaverage, outcomeaverage,other_punishmentaverage)
colnames(otherdata)[3] <- c("punishmentaverage")
#set the condition to -1 when it is other
otherdata$condition = -1

#create a variable that specifies the scenario
otherdata <- rownames_to_column(otherdata, var = "scenario")

#combine into a full dataset
selfother_averages.data <- rbind(selfdata,otherdata)

#mean center numerical variables
selfother_averages.data$actionavg_cent<-selfother_averages.data$actionaverage-(mean(selfother_averages.data$actionaverage))
selfother_averages.data$outcomeavg_cent<-selfother_averages.data$outcomeaverage-(mean(selfother_averages.data$outcomeaverage))
selfother_averages.data$punishmentavg_cent<-selfother_averages.data$punishmentaverage-(mean(selfother_averages.data$punishmentaverage))

#create interaction variables
selfother_averages.data$aa_x_cond <- as.numeric(selfother_averages.data$condition)*selfother_averages.data$actionavg_cent
selfother_averages.data$oa_x_cond <- as.numeric(selfother_averages.data$condition)*selfother_averages.data$outcomeavg_cent

#specify that categorical variables are factors
selfother_averages.data$scenario <- factor(selfother_averages.data$scenario)
selfother_averages.data$condition <- factor(selfother_averages.data$condition, labels = c("other", "self"))

```
#Self/Other descriptives
```{r descriptives}
punishment_desc <- ddply(selfother_averages.data, c("condition"), summarise,
               N    = length(punishmentaverage),
               mean = mean(punishmentaverage),
               sd   = sd(punishmentaverage),
               se   = sd / sqrt(N)
)
print(punishment_desc)

self_data<-selfother_averages.data[which(selfother_averages.data$condition=="self"),]
other_data<-selfother_averages.data[which(selfother_averages.data$condition=="other"),]

mosaic::cor(punishmentavg_cent~actionavg_cent,data=self_data)
mosaic::cor(punishmentavg_cent~actionavg_cent,data=other_data)

mosaic::cor(punishmentavg_cent~outcomeavg_cent,data=self_data)
mosaic::cor(punishmentavg_cent~outcomeavg_cent,data=other_data)
```

#Self/Other analyses
```{r linear models}

# full model
full_mod <- lme4::lmer(punishmentavg_cent ~ actionavg_cent + outcomeavg_cent + condition + aa_x_cond + oa_x_cond + (1|scenario),data=selfother_averages.data)

summary(full_mod)

#to test the hypothesis that the actionavg beta is significantly larger than the outcomeavg beta
library(car)
linearHypothesis(full_mod, c("actionavg_cent = outcomeavg_cent"))

#model without action main effect
NoAction_mod <- lme4::lmer(punishmentavg_cent ~ outcomeavg_cent + condition + aa_x_cond + oa_x_cond + (1|scenario),data=selfother_averages.data)

summary(NoAction_mod)
anova(full_mod,NoAction_mod)

#model without outcome main effect
NoOutcome_mod <- lme4::lmer(punishmentavg_cent ~ actionavg_cent + condition + aa_x_cond + oa_x_cond + (1|scenario),data=selfother_averages.data)

summary(NoOutcome_mod)
anova(full_mod,NoOutcome_mod)

#model without condition main effect
NoCondition_mod <- lme4::lmer(punishmentavg_cent ~ actionavg_cent + outcomeavg_cent + aa_x_cond + oa_x_cond + (1|scenario),data=selfother_averages.data)

summary(NoCondition_mod)
anova(full_mod,NoCondition_mod)

#model without action x condition interaction
NoActxCond_mod <- lme4::lmer(punishmentavg_cent ~ actionavg_cent + outcomeavg_cent + condition + oa_x_cond + (1|scenario),data=selfother_averages.data)

summary(NoActxCond_mod)
anova(full_mod,NoActxCond_mod)

#model without outcome x condition interaction
NoOutxCond_mod <- lme4::lmer(punishmentavg_cent ~ actionavg_cent + outcomeavg_cent + condition + aa_x_cond + (1|scenario),data=selfother_averages.data)

summary(NoOutxCond_mod)
anova(full_mod,NoOutxCond_mod)

```
#Self/Other plots
```{r new graphs}

# action vs punishment + outcome vs punishment divided by condition

# punishment vs action 
ggplot(selfother_averages.data,aes(y=punishmentaverage,x=actionaverage))+
  geom_point()+
  facet_wrap(~condition)+
  geom_smooth(method=lm)

# punishment vs outcome 
ggplot(selfother_averages.data,aes(y=punishmentaverage,x=outcomeaverage))+
  geom_point()+
  facet_wrap(~condition)+
  geom_smooth(method=lm)

#action vs. outcome
ggplot(selfother_averages.data,aes(y=actionaverage,x=outcomeaverage))+
  geom_point()+
  geom_smooth(method=lm)

#bar graphs

#punishment
ggplot(data=selfother_averages.data, aes(x=scenario, y=punishmentaverage, fill=condition)) +
  geom_bar(stat="identity", width=.9, position=position_dodge()) + scale_fill_brewer(palette="Set2")

#model without action or action x condition interaction
NoAct_NoActxCond.mod <- lme4::lmer(punishmentaverage ~ outcomeavg_cent + condition + oa_x_cond + (1|scenario),data=selfother_averages.data)
#get residuals from this model
selfother_averages.data$act_resids<-resid(NoAct_NoActxCond.mod)
#plot residuals predicted by action and condition
ggplot(selfother_averages.data,aes(y=act_resids,x=actionaverage))+
  geom_point()+
  facet_wrap(~condition)+
  geom_smooth(method=lm)

#model without outcome or outcome x condition interaction
NoOut_NoOutxCond.mod <- lme4::lmer(punishmentaverage ~ actionavg_cent + condition + aa_x_cond + (1|scenario),data=selfother_averages.data)
#get residuals from this model
selfother_averages.data$out_resids<-resid(NoOut_NoOutxCond.mod)
#plot residuals predicted by action and condition
ggplot(selfother_averages.data,aes(y=out_resids,x=outcomeaverage))+
  geom_point()+
  facet_wrap(~condition)+
  geom_smooth(method=lm)


```

```{r average study1 and study2 data}

#load the CSVs for study 1
study1_inout <- read_csv("~/Desktop/Studies/PunishmentAversion/Study2/study1_inout.csv",col_names = TRUE)

study1_selfother <- read_csv("~/Desktop/Studies/PunishmentAversion/Study2/study1_selfother.csv",col_names = TRUE)

study1_baseline <- read_csv("~/Desktop/Studies/PunishmentAversion/Study2/study1_baseline.csv",col_names = TRUE)

#pull apart data for averaging
study1_indata <- study1_inout[which(study1_inout$condition == "ingroup"), ]
study1_outdata <- study1_inout[which(study1_inout$condition == "outgroup"), ]
study2_indata <- inout_averages.data[which(inout_averages.data$condition == "ingroup"), ]
study2_outdata <- inout_averages.data[which(inout_averages.data$condition == "outgroup"), ]

study1_selfdata <- study1_selfother[which(study1_selfother$condition == "self"), ]
study1_otherdata <- study1_selfother[which(study1_selfother$condition == "other"), ]
study2_selfdata <- selfother_averages.data[which(selfother_averages.data$condition == "self"), ]
study2_otherdata <- selfother_averages.data[which(selfother_averages.data$condition == "other"), ]

#merge the study 1 and study 2 datasets into single datasets
full_baseline <- merge(study1_baseline, averages.data, by.x = "scenario", by.y = "scenario")
full_indata <- merge(study1_indata, study2_indata, by = "scenario")
full_outdata <- merge(study1_outdata, study2_outdata, by = "scenario")
full_selfdata <- merge(study1_selfdata, study2_selfdata, by = "scenario")
full_otherdata <- merge(study1_otherdata, study2_otherdata, by = "scenario")

#get study 1 and study 2 averages 
full_baseline$full_actionaverage <- rowMeans(full_baseline[ , c("actionaverage.x","actionaverage.y")], na.rm=TRUE)
full_baseline$full_outcomeaverage <- rowMeans(full_baseline[ , c("outcomeaverage.x","outcomeaverage.y")], na.rm=TRUE)
full_baseline$full_punishmentaverage <- rowMeans(full_baseline[ , c("punishmentaverage.x","punishmentaverage.y")], na.rm=TRUE)



full_indata$full_actionaverage <- rowMeans(full_indata[ , c("actionaverage.x","actionaverage.y")], na.rm=TRUE)
full_indata$full_outcomeaverage <- rowMeans(full_indata[ , c("outcomeaverage.x","outcomeaverage.y")], na.rm=TRUE)
full_indata$full_punishmentaverage <- rowMeans(full_indata[ , c("punishmentaverage.x","punishmentaverage.y")], na.rm=TRUE)

full_outdata$full_actionaverage <- rowMeans(full_outdata[ , c("actionaverage.x","actionaverage.y")], na.rm=TRUE)
full_outdata$full_outcomeaverage <- rowMeans(full_outdata[ , c("outcomeaverage.x","outcomeaverage.y")], na.rm=TRUE)
full_outdata$full_punishmentaverage <- rowMeans(full_outdata[ , c("punishmentaverage.x","punishmentaverage.y")], na.rm=TRUE)

full_selfdata$full_actionaverage <- rowMeans(full_selfdata[ , c("actionaverage.x","actionaverage.y")], na.rm=TRUE)
full_selfdata$full_outcomeaverage <- rowMeans(full_selfdata[ , c("outcomeaverage.x","outcomeaverage.y")], na.rm=TRUE)
full_selfdata$full_punishmentaverage <- rowMeans(full_selfdata[ , c("punishmentaverage.x","punishmentaverage.y")], na.rm=TRUE)

full_otherdata$full_actionaverage <- rowMeans(full_otherdata[ , c("actionaverage.x","actionaverage.y")], na.rm=TRUE)
full_otherdata$full_outcomeaverage <- rowMeans(full_otherdata[ , c("outcomeaverage.x","outcomeaverage.y")], na.rm=TRUE)
full_otherdata$full_punishmentaverage <- rowMeans(full_otherdata[ , c("punishmentaverage.x","punishmentaverage.y")], na.rm=TRUE)

#Recombine the data across conditions
full_inoutdata <- rbind(full_indata, full_outdata)
full_inoutdata$Condition <- factor(full_inoutdata$condition.x, labels = c("Ingroup", "Outgroup"))

full_selfotherdata <- rbind(full_selfdata, full_otherdata)
full_selfotherdata$Condition <- factor(full_selfotherdata$condition.x, labels = c("Other", "Self"))
```

```{r final plots}
ggplot(full_baseline, aes(x=full_actionaverage, y=full_punishmentaverage)) +
  geom_point(color = "black", alpha = .6, size = 3) +
  geom_smooth(method="lm", se = FALSE, color = "dark grey") +
  scale_y_continuous(limits = c(0, 100), breaks=c(0, 20, 40, 60, 80, 100), labels = str_wrap(c("Definitely Would Not Approve", "20", "40","60","80", "Definitely Would Approve"),10)) + 
  scale_x_continuous(limits = c(35, 60), breaks=c(35, 40, 45, 50, 55, 60), labels = str_wrap(c("Less Upset", "40","45", "50","55","More Upset"),6)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(size=16), axis.text.y = element_text(size=16), axis.title = element_text(size=22, vjust=0.35)) +
  xlab("How upset would you feel performing this action?") + 
  ylab(str_wrap("How likely would you be to approve of this punishment?", 30))

ggsave("Fig1a.jpg",height = 7, width= 11, dpi=600)

ggplot(full_baseline, aes(x=full_outcomeaverage, y=full_punishmentaverage)) +
  geom_point(color = "black", alpha = .6, size = 3) +
  geom_smooth(method="lm", se = FALSE, color = "dark grey") +
  scale_y_continuous(limits = c(0, 100), breaks=c(0, 20, 40, 60, 80, 100), labels = str_wrap(c("Definitely Would Not Approve", "20", "40","60","80", "Definitely Would Approve"),10)) + 
  scale_x_continuous(limits = c(30, 100), breaks=c(30, 40, 60, 80, 100), labels = str_wrap(c("Less Suffering", "40", "60","80", "More Suffering"),6)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(size=16), axis.text.y = element_text(size=16), axis.title = element_text(size=22, vjust=0.35)) +
  xlab("How much suffering would this action cause?") + 
  ylab(str_wrap("How likely would you be to approve of this punishment?", 30))

ggsave("Fig1b.jpg",height = 7, width= 11, dpi=600)

ggplot(full_inoutdata, aes(x=full_actionaverage, y=full_punishmentaverage, col=Condition)) +
  geom_point(alpha = .8, size = 3) +
  geom_smooth(method="lm", se = FALSE) +
  scale_y_continuous(limits = c(0, 100), breaks=c(0, 20, 40, 60, 80, 100), labels = str_wrap(c("Definitely Would Not Approve", "20", "40","60","80", "Definitely Would Approve"),10)) + 
  scale_x_continuous(limits = c(35, 60), breaks=c(35, 40, 45, 50, 55, 60), labels = str_wrap(c("Less Upset", "40","45", "50","55","More Upset"),6)) +
  theme_minimal() +
  theme(axis.text.x = element_text(size=16), 
        axis.text.y = element_text(size=16), 
        axis.title = element_text(size=22, vjust=0.35),
        legend.text=element_text(size=20),
        legend.title=element_text(size=20)) +
  scale_color_brewer(palette="Set2") +
  xlab("How upset would you feel performing this action?") + 
  ylab(str_wrap("How likely would you be to approve of this punishment?", 30))

ggsave("Fig2a.jpg",height = 7, width= 11, dpi=600)

ggplot(full_inoutdata, aes(x=full_outcomeaverage, y=full_punishmentaverage, col=Condition)) +
  geom_point(alpha = .8, size = 3) +
  geom_smooth(method="lm", se = FALSE) + 
  scale_y_continuous(limits = c(0, 100), breaks=c(0, 20, 40, 60, 80, 100), labels = str_wrap(c("Definitely Would Not Approve", "20", "40","60","80", "Definitely Would Approve"),10)) + 
  scale_x_continuous(limits = c(30, 100), breaks=c(30, 40, 60, 80, 100), labels = str_wrap(c("Less Suffering", "40", "60","80", "More Suffering"),6)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(size=16), 
        axis.text.y = element_text(size=16), 
        axis.title = element_text(size=22, vjust=0.35),
        legend.text=element_text(size=20),
        legend.title=element_text(size=20)) +
  scale_color_brewer(palette="Set2") +
  xlab("How much suffering would this action cause?") + 
  ylab(str_wrap("How likely would you be to approve of this punishment?", 30))

ggsave("Fig2b.jpg",height = 7, width= 11, dpi=600)

ggplot(full_selfotherdata, aes(x=full_actionaverage, y=full_punishmentaverage, col=Condition)) +
  geom_point(alpha = .8, size = 3) +
  geom_smooth(method="lm", se = FALSE) +
  scale_y_continuous(limits = c(0, 100), breaks=c(0, 20, 40, 60, 80, 100), labels = str_wrap(c("Definitely Would Not Approve", "20", "40","60","80", "Definitely Would Approve"),10)) + 
  scale_x_continuous(limits = c(35, 60), breaks=c(35, 40, 45, 50, 55, 60), labels = str_wrap(c("Less Upset", "40","45", "50","55","More Upset"),6)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(size=16), 
        axis.text.y = element_text(size=16), 
        axis.title = element_text(size=22, vjust=0.35),
        legend.text=element_text(size=20),
        legend.title=element_text(size=20)) +
  scale_color_brewer(palette="Set2") +
  xlab("How upset would you feel performing this action?") + 
  ylab(str_wrap("How likely would you be to approve of this punishment?", 30))

ggsave("Fig3a.jpg",height = 7, width= 11, dpi=600)

ggplot(full_selfotherdata, aes(x=full_outcomeaverage, y=full_punishmentaverage, col=Condition)) +
  geom_point(alpha = .8, size = 3) +
  geom_smooth(method="lm", se = FALSE) +
  scale_y_continuous(limits = c(0, 100), breaks=c(0, 20, 40, 60, 80, 100), labels = str_wrap(c("Definitely Would Not Approve", "20", "40","60","80", "Definitely Would Approve"),10)) + 
  scale_x_continuous(limits = c(30, 100), breaks=c(30, 40, 60, 80, 100), labels = str_wrap(c("Less Suffering", "40", "60","80", "More Suffering"),6)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(size=16), 
        axis.text.y = element_text(size=16), 
        axis.title = element_text(size=22, vjust=0.35),
        legend.text=element_text(size=20),
        legend.title=element_text(size=20)) +
  scale_color_brewer(palette="Set2") +
  xlab("How much suffering would this action cause?") + 
  ylab(str_wrap("How likely would you be to approve of this punishment?", 30))

ggsave("Fig3b.jpg",height = 7, width= 11, dpi=600)

```
